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ABSTRACT 

We study the evolution of embedded clusters. The equations of motion of the stars 
in the cluster are solved by direct N-body integration while taking the effects of stellar 
evolution and the hydrodynamics of the natal gas content into account. The gravity 
of the stars and the surrounding gas are coupled self consistently to allow the realistic 
dynamical evolution of the cluster. While the equations of motion are solved, a stellar 
evolution code keeps track of the changes in stellar mass, luminosity and radius. The 
gas liberated by the stellar winds and supernovae deposits mass and energy into the 
gas reservoir in which the cluster is embedded. We examine cluster models with 1000 
stars, but we varied the star formation efficiency (between 0.05-0.5), cluster radius 
(0.1-1.0 parsec), the degree of virial support of the initial population of stars (0-100%) 
and the strength of the feedback. We find that an initial star fraction M^^/Mtot > 0.05 
is necessary for cluster survival. Survival is more likely if gas is not blown out violently 
by a supernova and if the cluster has time to approach virial equilibrium during out- 
gassing. While the cluster is embedded, dynamical friction drives early and efficient 
mass segregation in the cluster. Stars of to > 2 Mq are preferentially retained, at the 
cost of the loss of less massive stars. We conclude that the degree of mass segregation in 
open clusters such as the Pleiades is not the result of secular evolution but a remnant 
of its embedded stage. 

Key words: star clusters ~ star formation - numerical simulations: hydrodynamic, 
N-body 



1 INTRODUCTION 

The loss of gas that signifies the end of star formation in 
an embedded pr oto-c l uster may also ca use the young clus- 
ter to dissolve l|Hillsl Il98(]| : iLada et all [l984i '). The energy 
output from all the massive stars typically exceeds the total 
binding energy of the embedded star cluster. The survivabil- 
ity of the cluster therefore depends on the efficiency with 
which the radiative, thermal and mechanical energy of the 
stellar outfiows couple to the inter-cluster gas ( and as long 
as th is energy is not carried away by outflows, I Dale et al.l 
l2005l ). The time scale on which gas is removed from the 
embedded cluster (and any spatial dependence of this pro- 
cess) affects the mass and number of stars in the surviv- 
ing cluster, but also its degree of mass segregation and the 
density profile. The parameters that strongly affect the out- 
come of the out-gassing phase are the star formation ef- 
ficiency (S FE, which by itself c o uld depend on feedback 
proce sses, iDale fc Bonneil l2008l : IVazguez -Semadeni et al.l 
I2OI0I ). the efficiency of the radiative coupling and the pop- 



* E-mail: polupes@strw.leidenuniv.nl 
© 2011 RAS 



ulation of the most massive stars in the proto cluster. 
Traditionally the critical value for the star formation ef- 
ficiency is take n to be SFE ~ 0.5 on th e basis of virial 
arguments (e.g. iHilld [lOSOl : lMathieulll983h . but later work 
showed that SFEs as low as 0.1 can still lead to bound clus- 



ters, albeit with a reduced mass jLada et al.l 198 4: Ad arnp 
'2000"; 'Gcvcr fc BurkertI I2001I : iBaumgardt fc Kroupal i2007[ 
I Goodwin, 2009l ) . O bserved star formation effi ciencies fall into 
the range 0.05-0.3 (| Williams fc McKed[l997l ) . It is therefore 
expected that clusters must experience considerable struc- 
tural disruption before settling. 

From an observational point of view, the number of 
embedded clusters is too high with respect to the num- 
ber of observed gas-free clusters. Within the solar neigh- 
borhood the observed distribution of cluster ages suggests 
that more than 95% o f embedded clusters dissolve into the 
field within 100 Myr (|Lada fc Ladall2003l ). The disruption 
of young clusters by the loss of the remnant gas provides a 
natural explanation, although alternatively the young clus- 
ters may be disrupted by tidal forces (Elm egreen fc Hunter! 
l2010l : lKruiissenll201ll ). In this scenario it is not the inter- 
nal evolution of the cluster that disrupts the structure, but 
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tidal shocks from passing gas clouds (see also iGieles et al.l 
|2006| ). In order to find out what the relative contribution is 
from these two processes we need to understand the effects 
of gas loss quantitatively. Ultimately the modeling of the 
formation of cluster populations is needed to interpret the 
observations as they play a pivotal role in galaxy evolution 
studies (e.g. ISchweizer &: SeitzeJ 1 19981 : iLarsen et al.l I2001I : 
iBastian et al.l 120051 ). Without understanding early cluster 
evolution we cannot relate the young cluster populations 
(which are easily observed to large distances) to the older 
stars and the stellar population in the field. Contrary to 
stars, clusters are sensitive to the galactic environment and 
if we would understand the starting conditions of the cluster 
population we could search for the signatures of later galac- 
tic evolutionary events such as galaxy mergers and major 
gas inflows. 

The expulsion of gas from young clusters has been 
studied by applying a time-dependent background po- 
tential within which the N-bod y syste m is integrated 
llLada et ai] Il984l: iGever fc Burk crt 2001 ; iBoilv fc Kroupal 
2003bl: iBaumgardt fc Kroupal 12007. : ,Chen fc Kol l2009l : 



Goodwinll2009l : ISmith et al.ll201ll '). The time scale for out 



gassing can then be controlled rather straightforwardly by 
the decay time scale of the background potential. This 
parametrization allows the study of different mass loss sce- 
narios withou t the complications of m odeling gas dynamics 
and feedback. iGever fc Burker^ (|200ll ) examined the evolu- 
tion of a cluster both using the external potential approxi- 
mation and for a limited set of combined N-body/ SPH simu- 
lations (using simplified feedback prescriptions). They found 
general agreement between the simulations with a dynamic 
gas component and a static gas potential, and found that 
relatively long timescales of gas loss are necessary for clus- 
ters to survive. To reconcile observed low SFEs with their 
simulations, they point out that although the global SFE 
may be low in a star forming region, locally the S FE may 
be m uch higher and allow for bound clusters (also lAdamd 
|2000| ). Furthermore they proposed sub-virial stell ar velocity 
dispe rsions as an alternative possibility (see also iGoodwinI 
l2009j ). 

IBaumgardt fc Kroupal (12007') conducted a grid of simu- 
lations in which they varied the SFE, gas expulsion timescale 
and external tidal field. They confirmed that clusters with 
a SFE as low as 10% could survive, although they found 
a strong dependence on the timescale (parametrized by an 
e-folding time) on which the gas potential was removed. Sur- 
viving clusters had expanded drastically; for SFEs of around 
25% this was a factor 3-4, nearly independent of the gas loss 
timescale. The velocity distributions became strongly radi- 
ally anisotropic in the cases where the gas was removed on 
timescales much shorter than a crossing time. 

It has become clear from observations and simulations 
that young stars start out in a structured , clum py, state 
l|Lada fc Ladal l2003l : iPortegies Zwart et all I2OI0I'). Recent 
work examining hiera rchical cluster formation ( Smith et al.l 
l201ll : lKruiissenf [ 20 111 ) suggests that sub-clustering increases 
the survivability. If clusters tend to be born sub-virial, they 
can survive even the loss of large fractions of intercluster gas, 
in whic h case the SFE may not be a g ood indicator of surviv- 
ability. [Smit^^£|al] (|201ll ) (see also iKruiissen,2011, ) found 
that the local gas fraction immediately prior to gas loss is 
a better diagnostic for the cluster survival than the global 



SFE. By studying sub-clusters in hydrodynamical simula- 
tions of colla psing mole c ular c louds where sink particles rep- 
resent stars. iKruiissenI (|201ll ) found that the environment 
around ~ 40 Mq clumps has a gas content of only 10%- 
30%. The relative lack of gas resulted from the accretion 
of gas onto stars (which depletes the gas reservoir) and the 
shrinkage of the virialized sub clumps (because the dynam- 
ics of these clumps decouples from the gas dynamics). For 
such low gas fractions no major disruption is expected when 
feedback clears the gas. On the larger s cales gas fractions 
are still considerably higher. lKruiisse"nl (|201ll ) argues that 
a similar effect will also cause the larger scales to become 
gas poor, if sufficient time elapses before the onset of feed- 
back. However, during this delay the different clumps would 
merge and relax, evolving from the highly structured state 
towards a smooth density distribution ( this w as also the 
case for the simulations for ISmith et al.l |20l3), where this 
happened within 3 Myr). The conditions at this stage would 
determine the outcome with regard to the final survival, and 
by this stage the local depletion of gas by accretion may no 
longer be present. 

The above studies have expanded considerably on the 
simple picture where clusters dissolution is determined solely 
or mainly by the single SFE parameter. Clusters seem to 
have a number of ways in which they can resist the dis- 
ruptive effect of the loss of their gas content: disruption is 
less or absent if the gas loss occurs on timescales longer 
than the crossing time or if the gas loss is delayed enough 
that shrinking due to relaxation processes can occur. On 
the other hand, if the stellar population of the cluster forms 
sub- virialized, some of the stars in the cluster can evolve 
into low energy orbits and these will not be ejected as easily, 
even if the SFE is low. Hence, to constrain the disruption 
process of clusters, we need to employ numerical simula- 
tions to unravel various complex physical processes at work. 
The determination of the gas loss timescale and the spa- 
tial distribution requires careful modeling of the response 
of the parent gas cloud to the mechanical luminosity of the 
young stars using realistic stellar evolution models (the feed- 
back from radiation, stellar winds and supernovae will be 
neither constant nor instantaneous) coupled to a hydrody- 
namic code. Following the evolution of the stellar compo- 
nent requires that the stellar dynamics be resolved in high 
detail, including two-body relaxation processes. For this an 
N-body code for coUisional gravitational dynamics is nec- 
essary. The relaxation processes also require that a real- 
istic initial mass function (IMF) is adopted, which is also 
required for the stellar evolution modeling. The stellar dy- 
namics and hydrodynamics code must be mutually coupled 
as they feel each other's gravity. Here we will present simu- 
lations that were conducted using the newly developed As- 
trophysical Multi-purpose Software Environment (AMUSE; 
IPortegies Zwart et al.ll2009l : |Pelupessv et al.ll201ll ) that cou- 
ple stellar dynamics of young clusters with the gas dynam- 
ics using realistic feedback physics. Using these simulations 
we will explore idealized but realistic scenarios of gas loss 
and analyze the trends with SFE, feedback strength and 
the dependence on the initial distribution. We will start by 
presenting the AMUSE package, the solvers and initial con- 
ditions of the runs we conducted in Section [S] The results 
are presented in Section [3] and discussed in|4] We summarize 
our main conclusions in Section (5) 



© 2011 RAS, MNRAS 000, - 



The evolution of embedded star clusters 3 



2 METHOD 

We study the evolution of embedded star clusters using a 
hybrid simulation environment. Our models combine the ef- 
fects of the self-gravity of the stars, their nuclear evolution 
and the hydrodynamics of the gas. The latter has contri- 
butions of the primordial gas content and of the gas lib- 
erated by the stars in stellar winds and supernovae. The 
various simulat ion ingredients are coupled self-consistently 
using AMUSE ijPortegies Zwart et al.ll2009l : |Pelupessv et al.l 
I2OIII ). In the code, a script in the AMUSE framework, we 
couple a coUisional N-body code, a Smoothed Particles Hy- 
drodynamics (SPH) code and a stellar evolution code to- 
gether to form a complete description of the physics that 
govern the evolution of an embedded star cluster. The inte- 
gration of the equations of motion of the stars was conducted 
using the fourth-or der predictor -corrector Hermite N-body 
solver phiGRAPE (jHarfst et al.l[2007, ). The dynam ics of the 
gas w as resolved using the SPH code Gadget-2 (jSpringell 
I2OO5I ). For the gravitational coupling between the stars and 
the self-g ravity of the gas we u sed the gravitational treecode 
Octgrav (jCaburov et al.ll201(]| ). The stars in our simulations 
were evolved usin g parametrized stel lar evolution tracks (us- 
ing the SSE code. lHurlev et aDbOOd ). Each of the numerical 
simulation codes is embedded in the AMUSE framework, 
and the code that realizes the self-consistent coupling be- 
tween each of the physical sub-domains turns out to be a 
relatively straightforward script. Within the AMUSE envi- 
ronment the individual sub-solvers are referred to as commu- 
nity codes. The algorithmes used by these community codes 
remain the same as described by their respective authors. 
The codes are written in different languages and have a wide 
range of programming styles. Because the adopted commu- 
nity codes are untouched, we will not further explain their 
operations, but simply refer the reader to the appropriate 
literature. The AMUSE framework that stitches each of the 
community codes together is one of the major ingredients 
that enables us to perform a consistent coupling between 
the various scales and physical domains. Therefore we give 
a short overview of the framework in the following section. 



2.1 AMUSE 

The AMUSE environment that we use to conduct the sim- 
ulations is a package that allows astrophysical codes from 
different domains to be combined to conduct numerical 
experiments (see Pelupessy et al. 2011 for a more com- 
plete description). AMUSE is a development of the Multi- 
physics and Multi-scale S oftware Environment (MUSE, 
IPortegies Zwart et al.ll2009l ) and is freely available for down- 
loacQ. The fundamental idea of AMUSE is the abstraction 
of the functionality of simulation codes into physically moti- 
vated interfaces that hide the complexity and numerical im- 
plementation of the codes. AMUSE presents the user with 
building blocks that can be combined into applications and 
numerical experiments. The binding language that stitches 
the codes together is Pythorfl, as the focus in these high 
level interactions is not so much performance (the computa- 
tional cost being concentrated in the component codes) but 



on algorithmic flexibility and ease of programming to allow 
rapid prototyping. 

An AMUSE application consists, roughly speaking, of 
a user script, one or more community modules and the com- 
munity code (fig[lj. The user script specifies the initial con- 
ditions and selects the simulation codes. The relationships 
between the community codes define the solver and for our 
case we construct a combined N-body/hydrodynamic solver 
using a direct N-body code and an SPH solver. The setup 
and communication with a community code is handled by 
the community module. This consists of an MPfl based com- 
munication interface with the code as well as unit handling 
facilities and an object oriented data model. The heart of an 
AMUSE application consists of th e solvers empl oyed: here 



these are the Gadget-2 SPH co de (ISpringelll2005t) . the Phi- 
GRAPE Hermite N-b odv code (iHarfst et al.ll2007h . the tree 
gravity code Octgrav (iGaburov et al.ll2010l ) and the stellar 



evolution code SSE (|Hurlev et al.ll2000l ) 



2.2 The combined solver 

An overview of the calling sequence of our combined solver 
is given in figure [51 It consists of a combined integrator for 
coupled gas/gravitational dynamics systems and a feedback 
prescription for mechanical energy input. 



2.2.1 BRID GE integrator 

The BRIDGE integrator l|Fuin et al.ll2007l ) provides a semi- 
symplectic mapping for the gravitational evolution in cases 
were the dynamics of a system can be split in different 
regimes. A typical application would be a dense star clus- 
ter in a galaxy where the internal dynamics of the former 
evolves on a relatively short timescale compared to the dy- 
namics of the ho s t gala xy. A similar idea was applied by 
ISaitoh fc Makinol (|2010l ) by splitting the gravitational and 
hydrodynamic evolution operators for simulating gas-rich 
galaxy mergers. They implemented the algorithm in a sin- 
gle monolithic code, whereas we adopt the concept of oper- 
ator splitting within AMUSE to couple different codes. In 
our case we will also employ a split between the gravita- 
tional evolution operator of the stellar component (which 
is governed by coUisional dynamics) and the hydrodynamic 
evolution operator (governed by SPH dynamics). 

The BRIDGE formulation is derived by applying 
the splitting method devel o ped in planetary dynam- 
ics (jWisdom fc Holmanl Il99ll : buncan et all 1 19981 ) to the 
Hamiltonian 



H 



= y y 



Gmiruj 
\\■r^-rj\\ 



(1) 



of a system of particles which consists of subsystems A and 
B by dividing the Hamiltonian into three parts 

H = Z^^+ 2^ - + 



, 2mi \\ri — rj 

Pi , GmirUj 



+ 



(2) 
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^ www . python . org 
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Figure 1. The AMUSE interface design. This diagram represents the way in which a community code ("code") is accessed from the 
AMUSE framework. The code has a thin layer of interfaces functions in its native language which communicate through an MPI message 
channel with the python host process. On the python side the user script ("AMUSE simulation script") only accesses generic calls 
("setup," "evolve" etc.) to a high level interface. This high level interface calls the low level interface functions, hiding details about units 
and the code implementation (together these form the community module of the code). The communication through the MPI channel 
does not interfere with the code's own parallelization. 



Gmiirij 



= Ha + Hb + H'Tb (3) 

where Ha and Hb are the Hamiltonians of subsystems A 
and B respectively and the cross terms are coUected in _ff'"'. 
The formal time evolution operator of the system can then 
be written as (H is the Hamiltonian vector field correspond- 
ing to a Ifamiltonian H): 

^tH _ ^T/2H'"'gT(H^+HB)gT/2H'"* 

The e^^™ operator represents momentum kicks because 
depends only on the position of the particles. The op- 
erator e'^''^-^^''^-' splits into secular evolution operators, be- 
cause Ha and Hb are independent. The evolution of the 
total system over a time step r can subsequently be calcu- 
lated from the mutual forces exerted by systems A on B and 
by kicking the momenta for time r/2, after which the two 
systems are evolved in isolation for a time r, after which 
the time step is finished by another mutual kick. This pro- 
cedure is beneficial if, for example, the time step allowed by 
the mutual interactions is considerably longer than that re- 
quired by the internal dynamics of the component systems, 
or if the gravitational dynamics of the two sub-systems need 
different solving strategies. 

Within AMUSE we view the various gravitational and 



hydrodynamical codes as implementations of the time evo- 
lution operators e'^"' (for this the interface specification pro- 
vides an evolve function). The momentum kicks are easy 
and reasonably fast to implement within the framework in 
Python once the forces are known, and for this the gravita- 
tional dynamics interface provides utility functions to query 
the gravitational forces at arbitrary positions. 

The gas dynamical SPH evolution can also be derived 
from a Hamiltonian formalism and thus the above directly 
carries over to a split between purely gravitational and 
SPH particles. The application of a BRIDGE type inte- 
grator to SPH is computationally efficient if the time step 
requirement for the hydrodynamics is more stringent than 
for the gravitational dy namics (which generally is the case, 
ISaitoh fc Makindl2010l ). Our application within AMUSE has 
the additional benefit that we can easily swap in and out dif- 
ferent gravitational solvers and SPH codes, even if they use 
completely different algorithms (such as in our case a direct 
Hermite code and a tree based SPH integrator). 

For validation, the resulting combined gravitational- 
hydrodynamical solver was checked in problems which could 
also be calculated using a monolithic SPH code. For this a 
number of static and dynamic test runs were conducted and 
the results where verified to be consistent with runs con- 
ducted purely with Gadget-2. First, in order to check that 
the alternation of hydrodynamic and gravitational evolution 
operators does not introduce artifacts in cases where grav- 
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Figure 2. The AMUSE gravitational/hydro/stellar evolution in- 
tegrator. This diagram shows the calling sequence of the differ- 
ent AMUSE elements in the combined gravitational/hydro/stellar 
solver during a time step of the combined solver. Circles indicate 
calls to the (optimized) component solvers, while rectangles indi- 
cate parts of the solver implemented in Python within AMUSE. 

itational and hydrodynamic forces are in equilibrium, we 
conducted static tests where the abihty of the code to main- 
tain a Plummer sphere in equiUbrium was examined. Tests 
where conducted using different particle numbers (10k-200k) 
and gas fractions (50 % and 100%). As a dynamic test we 
conducted the Evrard l|Evrardlll98"a ) collapse tests with the 
combined AMUSE solver. Both the dynamic and static tests 
showed satisfactory results compared to a conventional in- 
tegrated code jPelupessv et al.ll201ll ). 

2.2.2 Feedback prescription 

In our calculations we include feedback from stellar 
wind and supernovae. Radiative processes are not taken 
into account self-consistently, even though the ioniz- 
ing r adiation from massive stars can drive violent out- 
flows (|Rodriguez-Gaspar et al ]|l995l : lDiaz-Miller et al.lll998l : 
iDale et al.ll2005h : the relative strength of the radiative feed- 
back is expected to be of the same order of magnitude as th e 
mechanical feedback from stellar winds (|Fall et al.l l201(t ) , 
which is included in our calculations. 

Given these uncertainties (and the fact that we do not 
explicitly include radiative feedback) we have chosen to sim- 
plify the modeling by taking an adiabatic equation of state 
and not solving the thermal and chemical balance of the gas. 
We parametrize the energy budget in the absence of radia- 



tive losses with a coupling strength parameter which is 
the fraction of the total supernova and wind energy output 
that ends up as thermal energy in the ISM. The relative 
strength of the feedb ack is of the order of a few percent, the 
rest is radi ated away (' Wheeler et allll980l : ISiUch et al.|[l996l : 
iFrever et al. 2003 . 2006. ) . 

For the actual implementation of the feedback within 
our model we calculate the mechanical luminosity Lmech = 
Mw^ of the stars, from parametrized mass loss rates and 
terminal wind velocities (|Leitherer et al.lll992l : jPrinia et al.l 
[T99Q) . The mass liberated in the stellar wind of the stars is 
subsequently injected into the cluster medium. All gas par- 
ticles have the same mass rrigas throughout the simulation. 
The gas particles are released at rest relative to the star re- 
leasing them in a sphere with a radius equal to the smooth- 
ing length of the gas particles around the star, and the stellar 
mass is reduced accordingly. A consequence of this prescrip- 
tion is that each star loses mass in discrete steps of at least 
mgas. For the runs presented here the rrigas varies between 
3 X 10"^ - 6 X 10"^ M0. As a consequence the stellar wind 
feedback will exhibit some discreteness noise if the mass loss 
rate of a star is lower than ~ 7.5 x 10~* — 1.5 x 10~® Me/yr 
(Note that the feedback will still be approximately isotropic 
if a low number of particles is released due to the SPH 
smoothing). We verified that the results of our simulations 
with respect to the stellar response converges and therefore 
are not affected by this discretization. For the most active 
wind phases {M > lO~^M0/yr) the feedback prescription 
leads to 10s- 100s of particles released per timestep and a 
smooth wind (see section [5} . The internal energy u of the 
gas lost by the star is calculated from the total energy out- 
put of the star -Bmech since the last discrete mass loss event: 

u = /fb ^^22^ ^ I Lmcch / (rrrgrav - movo)dt(5) 

rrigrav rricvo "^^last 

here rrigi-av is the dynamical mass and rricvo the mass calcu- 
lated by the stellar evolution code (so reduced by the mass 
loss). In case that a star experiences a supernova explosion 
we inject Esn = 10^^ ergs in the escaping particles by added 
this amount of energy to iJmoch- The feedback efficiency pa- 
rameter /fb that accounts for the radiative losses is a free 
parameter, which we varied between /fb = 0.01 and 0.1. 
In order to track the rapid evolution caused by feedback 
we adopt a maximum time step for the gas particles of 10^ 
years and after supernova events the maximum hydrody- 
namic time step is reduced to 20 years for 30k years in or- 
der to preve nt time-stepping artefacts (a more elegant fix is 
descr ibed in lSaitoh fc Makindl2009l : lDurier fc Dalla Vecchial 
l201ll . but our approach has only a minor effect on the total 
running time) . We summarize the various numerical param- 
eters in Table [1] 

2.3 Initial conditions 

The simulated clusters in Table[2] are composed of a mix- 
ture of gas and stars (with A'^* = 1000) distributed in 
a Plummer sphere. Stellar-stellar gravitational interactions 
are smoothed (Plummer smoothing) with a smoothing 
length 0.001 x -Rscaic and the interaction between gas parti- 
cles and star particles are smoothed (using the SPH spline 
kernel) with the gas particle smoothing length. The stel- 
lar masses are assigned using a Salpeter IMF between 0.1 
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Table 1. Overview of the numerical parameters: Atgas is the maximum time step for hydrodynamics, Atf^ 
defines the interval between feedback, Atbridge is the coupling time for mutual gas - star gravity kicks, tgim is 
the total simulated time, Csoft is the gravitational smoothing adopted (for star-star interactions. Interactions 
with gas particle s are smoothed using the local smoothing lengt h) , 77 is the usual N-body time step parameter 
for PhiGRAPE l lAarsethll20oi) and the lSarnes fc HulJ l IlQSd) opening angle for Octgrav. 



parameter; 


Aigas 


Aifb 


^^bridgc 




value: 


0.001 Myr 


0.04 Myr 


0.04 Myr 


30 Myr 0.001 i?,scalc 0.01 0.5 



Table 2. Overview of initial conditions. 



Model 




SEE 




Afgas 


Afgas 




/fb 


comment 








(M0) 




(Mq) 


(pc) 






Al 


1000 


0.05 


355 


100k 


6740 


0.5 


0.1 




A2 


1000 


0.3 


350-366 


100k 


815-855 


0.5 


0.1 


x4 


A3 


1000 


0.5 


355 


100k 


355 


0.5 


0.1 




A4 


1000 


0.05 


355 


100k 


6740 


0.5 


0.01 




A5 


1000 


0.3 


355-366 


100k 


815-855 


0.5 


0.01 


x4 


A6 


1000 


0.5 


355 


100k 


355 


0.5 


0.01 




A7 


1000 


0.3 


355 


100k 


840 


0.2 


0.1 




A8 


1000 


0.3 


355 


100k 


840 


0.2 


0.01 




A9 


1000 


0.3 


355 


100k 


840 


1.0 


0.1 




AlO 


1000 


0.3 


355 


100k 


840 


1.0 


0.01 




Bl 


1000 


0.05 


355 


100k 


6740 


0.5 


0.1 


sub-virialized 


B2 


1000 


0.05 


355 


100k 


6740 


0.5 


0.01 


sub-virialized 


B3 


1000 


0.05 


355 


100k 


6740 


1. 


0.1 


sub-virialized 


B4 


1000 


0.05 


355 


100k 


6740 


1. 


0.01 


sub-virialized 



and 100 M0, with the additional constraint that the most 
meissive star in each cluster has a mass close to the me- 
dian majcimum mass for a cluster of this size (« 22 M©). 
This additional requirement is co nsistent with the expected 
maximum mass for small clusters ijWeidner &: Kroupall200^ : 
IPortegies Zwart~et al. 2010). The masses of the stars are as- 
signed independently of the position in the cluster (constant 
SFE), so no initial mass segregation is imposed. Stellar col- 
hsions are not taken into account and during the simulation 
no gas accretion unto existing stars or new star formation is 
allowed. Note that any accretion of gas onto the stars would 
be very limited since we start from a smooth gas distribution 
and no cooling is included. Feedback induced star formation 
can also not occur in our approac h, but is not very l ikely in 
the systems we simulate (see also lDale fc Bonnelj|2008l ). 

In the series of calculations indicated with the letter "A" 
in Table[2] we explore the effect of the formation efficiency, 
the initial scale length and the feedback efficiency. We run 
models with a SFE of 0.05, 0.3 and 0.5. The initial Plummer 
radius of the initial conditions is either 0.1, 0.5 or 1.0 pc. 
Simulations with a feedback efficiency of /fb = 0.01 and 
0.1 are performed. For small clusters the number of high 
mass stars can vary quite substantially between different 
realizations of the IMF, and we perform the model A2 and 
A5 with realizations of the IMF generated from different 
random seeds to examine this effect. We have performed 
additional simulations in which we varied the number of gas 
particles to tests that the results are independent of the 
resolution of the gas dynamics. 

We also run a series of models in which the stars are dy- 
namically cold (these models with a sub-virialized velocity 
dispersion are designated as model set "B" in Table [2|. This 



choice is motivated by the fact that the initial stellar veloc- 
ity dis persion comes from tur bulent velocities in the parent 
cloud (|Gever fc BurkertHiooH ). The stellar component may 
form with a lower velocity dispersion than that required to 
be dynamically stable if the turbulence in the parent gas 
cloud has a scale length smaller than the Jeans length or if 
the gas cloud is (partially) pr essure supported o r partially 
supported by magnetic fields (| VerschuerenI 1 1990l ) . For this 
set of simulations we decrease the stellar velocities to 20% 
of the virialized values. 



3 RESULTS 

In Figure [3] we show the stellar and gas distribution of our 
'baseline' run A2 and its low feedback equivalent, run A5. 
At four moments in time we plot the stars and a cut through 
the 2 = plane of the gas density distribution. Looking at 
the first A2 frame (at 0.96 Myr) we see the early stages 
where stellar winds blow bouyant bubbles that rise in the 
potential of the star cluster. As the mechanical luminosity 
increases these bubbles grow until they start blowing away 
sizable fractions of the cluster medium and a free-fiowing 
wind develops (4.37 Myr frame). The strong feedback pro- 
ceeds to unbind most of the gas of the cluster. At approxi- 
mately 9.5 Myr the cluster ISM has been blown away - the 
gas that is visible in this frame originates from the strong 
AGB wind of the heaviest progenitor (21 Mq). The super- 
nova that explodes at 9.54 Myr leaves the cluster devoid of 
gas. 

For the A5 simulation (/fi^ = 0.01) (both simulations 
use the same initial realizations for the IMF and stellar po- 
sitions) the initial wind stages proceed less violently, with 
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Figure 3. Stellar and gas distribution of the A2 and A5 run. Left panels show the gas and stellar distribution of the A2 run, Right panels 
of the A5 run. Snapshots are labelled with their time in the lower right corner. Shown as a density plot is a slice through the midplane of 
the gas density. Points show the stars in 4 mass groups: m* 0.9 (smallest red dots), 0.9 Mq ^ 2.5 Mq (intermediate yellow 

dots), 2.5 Mq ^ rrii, ^ 10 Mq (intermediate green dots) and ^ 10 Mq (large light blue dots). 
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Figure 4. Lagrangian radii and core radius for runs Al-AlO. Plotted in red dashed lines are the time evolution of the 
(0.02,0.05,0.1,0.2,0.5,0.75 and 0.9 ) Lagrangian mass radii of the stellar mass distribution. Blue line is the core radius. Vertical dot- 
ted lines indicate supernova events. 



smaller bubbles (0.96 Myr), while a free-flowing wind does 
not develop (compare the frame at 4.37 Myr) until just be- 
fore the supernova. The main difference between the A2 and 
A5 run is that most of the cluster gas remains present in the 
latter case, until the first supernova, at which stage the re- 
maining gas is blown away. Just before the supernova the 
cluster is much more compact than in the A2 run. 

The corresponding evolution of the Lagrangian radii 
and core radius and central gas density (of models Al-AlO) 
is shown in Figure [3] and we plot in Figure [S] the (instanta- 
neous) bound fraction of the gas and star component. For 
run A2, in the early phase (before 8 Myr) the gas is lost by 



stellar winds. From the beginning of the simulation the gas 
starts to boil off vigorously under the influence of the stellar 
winds. The stellar distribution reacts by expanding. Within 
3 Myr the core radius increases a factor « 2.5 (within this 
time approx 2/3 of the gas is lost) and stays more or less 
constant thereafter. This is also visible in the maps of Fig.|31 
where after the 4.37 Myr frame the stellar distribution re- 
mains qualitatively similar. Fig. [5] also shows that when the 
first supernova goes off, the gas in the cluster has already 
been completely removed. Later supernovae affect the stel- 
lar distribution very little. If we look at the bound fraction 
of stars shown in Figure [S] we see that « 60% of the stars 
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Figure 5. Stellar and gaseous bound mass fractions for runs Al-AlO. Plotted is the bound stelllar fraction (blue line) and the bound 
gaseous fraction (red dashed line) as fraction of the total initial system mass. Vertical lines indicate supernova events. 



remain bound at the end of the simulation. The escape of 
stars happens mostly at the time of early gas loss, but with a 
noticeable delay. Note that the bound fraction plotted here 
and for the other runs is an upper limit because this fraction 
is determined from the instantaneously bound stars (so at 
a particular point in time) and does not account for further 
unbinding by the escape of marginally unbound mass (this 
can account for part of the observed delay). 

Comparing runs with different SFE (models A1,A2 and 
A3) we see that in the case of low SFE (run Al) most of the 
gas is retained up until the moment of the first supernova 
(see fig. [5|. The gas completely unbinds because of the su- 



pernova at 9.54 Myr, and the cluster is essentially dissolved 
immediately at that point. The high SFE case A3 is very 
similar to the A2 case. The initial gas loss may be some- 
what steeper, but in this case (fig. [S]) the loss of stars of the 
cluster proceeds more slowly. Looking at the evolution be- 
fore the first SN, it is clear that the greater gas mass (20 x 
more) gives the Al run greater resilience to the disturbance 
from stellar winds. Considering the escape of stars this leads 
to counteracting effects: when the fraction of gas is large, the 
loss of gas can more easily disturb the stellar distribution. 
However this loss happens initially much more slowly. In 
the A2 and A3 runs on the other hand, the gas starts to be 
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blown away right from the beginning, but the fraction of gas 
is small enough not to cause major disruption for the stel- 
lar distribution. When the supernova explodes, the gas of 
the Al run is instantaneously removed, leading to the quick 
dissolution of most of the cluster. On the other hand, about 
70% of the stars is retained in the A3 run. 

If we look at the runs with the lower effective feedback 
strength (A4, A5 and A6), we see that in the SFE = 0.05 
case (A4) gas is still present after the first and second su- 
pernova. The stellar component reacts completely differently 
for the A4 run (compared with Al): in this case 50% of the 
stars remain bound initially, where only 10% remains in the 
Al run. In the plot of the Lagrangian and core radii of run 
A4 (fig. 3)) the core radius initially increases after the super- 
nova, but then grows smaller again, as the cluster relaxes - 
markedly different from the continuing expansion in the Al 
case. For the other two runs (A5 and A6), the lower effective 
feedback strength mainly results in a much less violent boil 
off compared with the high /fb cases (e.g. A2 vs A5) . Gas is 
still present at the time of the first supernova. For the stel- 
lar distribution this means that initially the stellar content 
remains intact (in particular the A5 run shows an almost 
constant bound fraction for the stellar component up until 
9.5 Myr). The final bound stellar fractions and Lagrangian 
radii are qualitatively similar. The reason for this is that in 
both cases the total feedback energy is sufficiently large to 
unbind the cluster gas eventually. Much of the extra energy 
in the simulations with high feedback is carried away by the 
escaping gas. In the high feedback case gas loss by stellar 
winds is much more important, whereas in the low feedback 
case gas is driven out by supernovae. The final stellar distri- 
bution of the final cluster is little affected by the feedback 
efficiency over the range explored here, unless the gas frac- 
tion is very high. 



3.1 Different realizations 

For the low mass clusters we examine here different real- 
izations of the stellar density distribution and the IMF will 
have very different masses and orbits for the heavy stars (and 
thus in the number, timing and location of supernovae) just 
due to low number statistics (the number of stars larger than 
10 Mq is ~ 4). We have run additional realizations of the A2 
and A5 simulations to quantify this effect. For the A2 runs 
the range in final stellar bound fractions is 50 — 70%, with 
the final core radii lying between 0.6 and 1.2 parsec. For 
the A5 runs we find similar ranges, 50 — 65% and 0.8 — 1.5 
parsec. However if we compare the runs in more detail we 
see considerable differences. Two examples of this are given 
in figures |6] (Lagrangian and core radii) and [7] (bound frac- 
tions). The A2c run quickly loses 60% of its mass within 
1.5 Myr, after which the gas loss slows down momentarily. 
This results in a stellar distribution that quickly expands, 
but then contracts again until the remaining gas is expelled. 
The A5d run seems to follows a similar pattern to the A5 
run, if one compares the core radius. However if one looks 
more closely to the inner Lagrangian radii the A5d run is 
very different from the A5 run. In the latter case the (e.g.) 
10% Lagrangian radius hovers around 0.2 parsec after 10 
Myr, while in the A5d run it expands more than 4x further, 
to ~ 0.9 parsec. 
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Figure 6. Lagrangian radii, core radius and central gas density 
for the alternative runs A2c and A5b. These differ from the corre- 
sponding runs in Fig.|3]only in the random seed that was used to 
sample the density distribution and IMF. Plotted in red dashed 
Unes are the (0.02, 0.05, 0.1, 0.2, 0.5, 0.75 and 0.9) Lagrangian 
mass radii of the stellar mass distribution. Blue line is the core 
radius. Vertical dotted lines indicate supernova events. 
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Figure 7. Stellar and gaseous bound mass fractions for the al- 
ternative runs A2c and A5b. These differ from the corresponding 
runs in Fig. |4] only in the random seed that was used to sample 
the density distribution and IMF. Plotted is the bound stelUar 
fraction (blue line) and the bound gaseous fraction (red dashed 
lines) as fraction of the total initial system mass. Vertical lines 
indicate supernova events. 
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Figure 8. Lagrangian radii and core radius for runs B1-B4. Plotted in red dashed lines are the (0.02, 0.05, 0.1, 0.2, 0.5, 0.75 and 0.9 ) 
Lagrangian mass radii of the stellar mass distribution. Blue line is the core radius. Vertical dotted lines indicate supernova events. 





Figure 9. Stellar and gaseous bound mass fractions for runs B1-B4. Plotted is the bound stelUar fraction (blue lines) and the bound 
gaseous fraction (red dashed lines) as fraction of the total initial system mass. Vertical lines indicate supernova events. 



3.2 Sub-virial velocity dispersions 

The time evolution of the Lagrangian radii and the bound 
fractions of the models with sub-virialized velocity disper- 
sion are plotted in figures [8] and |9] Figure [8] shows that the 
stellar distributions quickly shrink compared to the viri- 
alized cases: e.g. for Bl the 0.5 Lagrangian radius drops 
roughly a factor 2 within less than 0.5 Myr. This collapse 
of the stellar distribution can increase the survival of the 



stellar cluster, bec ause the collapse increases the local SFE 
l|Smith et al.ll20lll ). The Bl model shows this effect, because 
compared with the Al model it shows a much less drastic 
expansion of the core radius (fig. [S]), although in the end 
the cluster still becomes completely unbound (fig. [9)l . How- 
ever, surprisingly, if we look at the B2 model with a lower 
/tb (which is the sub-virialized version of the A4 run) we 
see exactly the opposite: in this case the cluster dissolves 
quickly, where as one may expect on the basis of the Al and 
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Figure 10. Stellar mass fraction at the (stellar) half mass radius. 
Plotted is the mass fraction at the half mass radius of the stellar 
distribution as a function of time (in the initial stages of the 
simulations) for the sub-virial velocity dispersions, runs Bl (red) 
and B2 (green), and their equivalent with full velocity dispersion 
Al and A2 (blue and cyan dotted lines). 



A4 run that the B2 run would dissolve much more slowly. 
Runs B3 and B4 show a similar pattern. The reason for this 
becomes clear if we examine figure [TOl where we have plotted 
the stellar mass fraction (the ratio of stellar mass to total 
mass) within the stellar half mass radius as function of time 
up to the time of the first supernova. The dotted lines here 
show the Al and A4 runs, while the drawn lines are the Bl 
and B2 runs. The Bl and B2 runs initially show an increase 
of a factor of 2 in stellar fraction (which is smaller than ex- 
pected from the contraction of the half mass radius shown in 
fig- m indicating the gas distribution also contracts). After 
that they both show a slow increase, but crucially the Bl 
run shows a stronger increase to around 23% vs the B2 run. 
Although the increase is modest, it is probably the cause of 
the differ ence between the Bl and B2 r un. This is consistent 
with the Baumga rdt fc Kroupal (|2007l ) results, which show 
a very sensitive dependence of the bound fraction on SFE 
around a value of 0.25. 

Note that Figure fTOl shows that the stellar fraction of the 
A2 run is almost the initial value (SFE= 0.05). In this case 
the survival of the cluster after the first supernova is not due 
to a high stellar mass fraction, but due to the fact that not 
all gas is expelled (fig. [5]). The effective SFE hence increases 
to SFE « 0.05 + 0.2 (measured after 1 Myr) if one takes 
into account the gas that is not expelled. Hence the Bl and 
A2 run show an increased survivability due to completely 
different mechanisms: the Bl stellar component clears its 
environment of gas slowly through stellar winds and adapts 
to the slowly changing potential, while the stellar component 
of the A2 run is held together after the first supernova by 
remnant gas (in spite of this only being a very minor fraction 
of the initial gas component). 
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Figure 11. Mass fraction versus number fraction for the A2 run. 
Plotted is the mass fraction at given radii versus vs the num- 
ber fraction at those same radii (at 30 Myr). A cluster without 
ma ss s egregation would follow the dotted diagonal line. FoUow- 
ing lCo nvcrsc & Stahlcr ( 2008) we take the area between the curve 
and the diagonal as a measure for the mass segregation (the Gini 
coefficient for the cluster). 
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Figure 12. Plot of the Gini coefficient G as a function of time 
for different runs. The solid lines (red, green and blue) are runs 
with different initial half mass radius (A5, A8 and AID runs re- 
spectively). Red dashed and dotted lines are the A2 run and A6 
run. G is the area (x2) between the graph of mass fraction vs 
number fraction and the diagonal (Fig. Illl l. 
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Figure 13. Cumulative bound number fraction as a function of mass. Plotted is, as a function of mass, the number fraction of stars 
heavier than that mass that are still bound at 10, 20 and 30 Myr. Shown are models from the basis set Al-AlO which result in a bound 
cluster (and the A2c run). Note that errors on these are not plotted - these increase from left to right from 0.03 to 1 for the right most 
point (because the right most point is based on a single data point). 



3.3 Mass segregation and selective mass loss 

Visual inspection of the maps such as those in figure [3] shows 
that the resulting remnant cluster are mass segregated. A 
useful way to quanti f y the m ass segregation was presented by 
[Converse &: Stahleil ()2008l ) - they plotted the number frac- 
tion enclosed at given radii /jv against the enclosed mass 
fraction at those radii (/m). For mass segregated clusters 
(e.g. the Pleiades cluster in their paper) /a/ will lie above 
/jv (an example of this procedure for our simulation A2 is 
given in fig.lll[l. A measure of the mass segregation can then 
be given by calculating the area between the graph of /m 
and line of equality, which Converse & Stabler labeled the 
Gini coefficient G for the cluster. 

We have plotted G for various runs as a function of time 
in figure [121 Looking at run A5 (solid red line) we see that 
in the beginning G rises quite quickly, within 5 Myr; after 
that G is more or less constant, G ~ 0.2 — 0.3. At the end 
of the simulations the cluster is hence in a mass segregated 
state. A similar picture emerges for run A6 (higher SFE, 
dotted line in Fig. I12|) and for the higher feedback run (A2, 
dashed line). In the latter case G plateaus at a somewhat 
lower value. The other two drawn lines in figure[T2]show the 



runs with a higher and lower initial half mass radius, runs 
AlO and A8. As is expected in these cases the initial rise in 
G is slower (AlO) and faster (A8) respectively. In all cases 
the final G settles around G ft^ 0.2. 

It is interesting to note that [Converse fc" StahleJ (|2010l ) 
tried to model the Pleiades cluster in some detail by run- 
ning pure N-body models, and they could not explain the ob- 
served G unless they started their simulations with a nonzero 
Gini coefficient of at least G ~ 0.2. Although our models do 
not try to specifically reproduc e the Pleiades cluster, thi s 
seems consistent with Figure [T2] [Converse fc StahleJ (|2010l ) 
started their simulation after the dispersion of gas and our 
models suggest that the nonzero G is a remnant of the fast 
evolution before gas dispersal. After gas dispersal the clus- 
ter is left in a much mor e dispersed state (the models of 
[Converse fc Stahleil (|2010l ) had virial radii ^ 3 pc), with a 
G close to its present value. This suggests that the mass seg- 
regation in open clusters such as the Pleiades is not a result 
of the secular evolution of its present configuration, but a 
remnant of its embedded evolution stage. 

The fact that mass segregation in our simulations oc- 
curs early in the embedded stage of the cluster evolution also 
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has a consequence for the stellar population of the remnant 
cluster. In Figure [13] we plot the number of stars heavier 
than a given mass which end up in the remnant cluster (i.e. 
are bound) as a function of star mass. The panels in fig- 
ure [13] show that in most cases where there results a bound 
cluster, the heavier stars are preferentially retained (11 out 
of 15 runs show this pattern). A number of models do not 
show this pattern, especially A7 and A8. Note these clus- 
ters have drastically expanded. Note also that while the A7 
and A8 runs do not show strong preferential retention of 
heavy stars, they do show mass seg regation. It was previ- 
ously noted (jBoilv fc Kroupall2003al ) that simple virial the- 
orem arguments underestimate the size of clusters after gas 
clearing, because there are always a number of stars in low 
velocity orbits. The preferential retention of heavy stars can 
be understood by noting that the mass segregation during 
the initial evolutionary stages happens because of the par- 
titioning of stars in equal energy orbits: the heavy stars end 
up in low velocity orbits that are much more difficult to 
unbind. The heavy stars that are more easily retained may 
skew the IMFs that are derived from intermediate age clus- 
ters (as long as they are not observed in the initial stage) 
and this may result in an observable signature. Note that 
this could be tested for the Pleiades cluster, which is al- 
ready known to have an e xcess of heavy stars in the center 
ijConverse fc StahleillioTol . e.g.). 



3.4 Ratio of age and crossing time 

Loosely bound associations can be distinguished from 
'true' stellar clusters using limited observables by consid- 
ering the ratio 11 of age Td a nd the crossing time Tcross 
l|Gieles fc Portegies Z'TOrtll201ll 'l: 

n = Tcl/Tcross (6) 

where the Tcross is derived from 

- ^° (cf ) 

where the effective radius R^s and mass AI are easily esti- 
mated observable quantities. We determined the time evo- 
lution of n using Eq. [7] for our simulations by measuring 
the iioff and M. The resulting 11 are plotted in figure 1141 
Note that we plot the 11 for the case where we identify the 
cluster as being composed of all stars initially present or 
alternatively where only the bound stars are taken into ac- 
count. The evolution of the 11 of the clusters in Figure [14] 
can be understood as follows: in the beginning 11 increases 
linearly but as the cluster is disturbed the IT ratio is "reset," 
as the Res increases and the mass decreases. The tracks of 
our clusters end mostly with either a 11 « 1 (bound case) 
or with a very low 11 in which case they should be dissolved 
(and no longer recognized as clusters). This is consistent 
with a possible explanatio n for the distribution of II as a 
function of mass noted by l|Gieles fc Portegies Zwar tl l201ll ). 
They found that for young stellar associations both low and 
high n values are present. While older stellar associations 
were exclusively found with 11 ^ 1 (and identified as proper 
clusters) . 



4 DISCUSSION 

We performed simulations of relatively small clusters of 1000 
stars with gas. Our simulations self-consistently take the dy- 
namical evolution of the stars, their internal nuclear evolu- 
tion and the dynamics of the embedded gas and the gas liber- 
ated by stellar evolution into account. Each of the physical 
domains, stellar dynamics, stellar evolution and hydrody- 
namics, is resolved with well tested numerical solvers which 
have been developed independently. The coupling between 
the domain-specific solvers is realized with the AMUSE 
framework, which enables us to couple codes written in a 
wide variety of languages in a homogeneous, transparent and 
self-consistent fashion. The final code has the form of a rel- 
atively simple Python script, which is general in its purpose 
of simulating embedded star clusters. 

Our model clusters start with a mix of stars and gas in 
virial equilibrium at a hypothetical moment we call t — 0, 
at which all stars are on the zero-age main-sequence. The 
distribution of gas and stars was taken identical with a star- 
formation efficiency to control the relative mass fraction of 
the gas with respect to the stars. The feedback strength is 
parametrized using a parameter /fb. We constrain the ini- 
tial mass function to follow a Salpeter IMF with a minimum 
mass of 0.1 Mq and an upper limit of 100 Mq, but we fixed 
the maximum mass of a cluster star to 22 Mq in order to 
prevent random fluctuations from the coincidental genera- 
tion of a more massive star, because it turns out that the 
most massive star in our cluster dominates the out-gassing 
by its efficient segregation to the cluster center, the strong 
stellar wind and the moment the star experiences a super- 
nova. 

In our simulations a star formation efficiency SFE ^ 
0.05 results in a boun d cluster, which is consisten t with 
earlier simulations (e.g. iBaumgardt fc Krou"pall2007l ). How- 
ever, a more detailed analysis of our simulation results re- 
veals that the process is considerably more complicated. It 
turns out that the qualitative result with respect to final 
bound fraction can be approximated reasonably well with 
the parametriz ation of the SFE, but that th e entire process 
is much richer. iBa um gardt fc K roupal (200'f) found that the 
most important parameters affecting the final cluster were 
the SFE and the gas loss timescale. Our simulations con- 
vincingly demonstrate that the simple parametrization in 
terms of the SFE and/or gas loss timescale is a gross over- 
simplification of the process of stellar-supported out-gassing 
of the young cluster. For example for the low feedback A4 
run the initial supernova is not enough to disperse all cluster 
gas. After the first supernova 15% — 20% of the gas remains 
bound. From the standpoint of the dynamics of the stars this 
remnant gas increases the effective SFE. Hence the cluster 
initially relaxes (between 10-20 Myr) with a « 50% bound 
fraction. Comparing with the Baumgardt fc Kroupa results 
(Figure 1), we see that this is a reasonable value if we con- 
sider their SFE= 0.2 curves: at these SFE a bound fraction 
of 50% is quite possible (even not accounting for any radial 
dependence of the gas loss). Further supernova events blow 
away the remaining gas, but do so only after the cluster has 
adapted itself dynamically. On the other hand, the Bl run 
shows if the gas loss occurs more gradually, the cluster can 
also adapt to gas loss (albeit in our simulations not enough 
to give rise to a bound cluster). 
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Figure 14. Plot of 11 = age/Tcr for models A1-A6. Tcr is determined as in Eq. where Tcr depends on i^^ff s,nd M 
jGieles Portegies Zwart||201l| ). Two different cases are plotted: taking into account bound stars only (red line) or alternatively in- 
cluding all the stars (blue dashed lines). 



The time scale on which the gas is expelled cannot be 
described as a simple function, but depends on the distri- 
bution of the masses and positions of the stars, when they 
start blowing a strong stellar wind and the details regarding 
the moment of the supernovae (and the adopted feedback 
strength). Considering the wealth of variation in cluster pa- 
rameters that result from this complex process the question 
remains why young star clusters show a relatively narrow 
range of properties, whereas we would naively expect a much 
wider range of cluster radii, and masses. At this moment we 
are unaware what key parameter may play a moderating 
role in the range of observables for young star clusters. Note 
that we did not take into account external limitations, such 
as the galactic tidal fields or the environment of the parent 
cloud. 

In our simulations we do not model the stellar feed- 
back fully self-consistently. The main uncertainties in our 
prescription concern the uncertainty in the efficiency of the 
coupling of the mechanical luminosity, the constraints of the 
numerical method and the lack of radiative feedback in our 
models. In future simulations we expect to deal with this 
latter point by resolving the radiative feedback using a ra- 
diative transfer code, but for now this remains part of the 
uncertainty encoded in the free parameter for the feedback 
strength we introduced in our models. This efficiency with 
which the output of stars couples to the sur rounding gas 
may also d epend on the mass of the star itself (jFrever et al.l 
I2OO3I . [ 20061 ). and on the state of the gas: the density, tem- 
perature but also on its dumpiness (an aspect we have not 
modeled). We do demonstrate, however, that the core radii 
and mass of the surviving cluster does not depend strongly 
on this efficiency parameter. 

On the numerical side, it is well known that SPH 
shows relatively poor performance in the regime of explo- 



sive shock waves: for example the peak densities are gener- 
ally smoothed out. This could be important for the energy 
budget of the feedback, as the cooling happens mainly in 
these dense regions. On the other hand, the speed and post- 
shock states are well reproduced (|Springelll2005l : iPelupessvl 
l2005h in SPH, so with the limitations we imposed on the 
details of the thermal evolution (the fact that we did not 
try to account self-consistently for the energy budget of the 
feedback) this has probably not a major impact. In future 
work, AMR employing Godunov-type solvers for the gas dy- 
namics may be employed to improve this aspect. This would 
necessitate very high resolution simulations to reach a con- 
verged solution and the inclusion of much more physics to 
follow the (non-equilibrium) heating and cooling processes 
in the gas. 



5 CONCLUSIONS 

• We find that a SFE > 0.05 is necessary to result in 
bound clusters. 

• The actual dispersion of clusters at SFE= 0.05 may oc- 
cur gradual or sudden depending on the details of feedback. 
Before eventual dissolution, stable phases in the cluster evo- 
lution may occur if not all gas is expelled. Contributing to 
this is the increase in the effective SFE from the point of 
view of the stellar dynamics if feedback does not blow away 
all gas. 

• We find that statistical variations can have a consider- 
able effect on the internal structure of the remnant cluster. 

• We find different modes of gas loss depending on 
the feedback coupling strength adopted. A high coupling 
strength results in a gradual loss of gas on timescales of « 5 
Myr. For the low feedback coupling stellar winds are not 
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able to disperse the gas: the timescale of gas loss in this case 
is very short (due to supernova explosions). 

• More detailed modeling is necessary to put constraints 
on the energy budget and thermal evolution of the gas, espe- 
cially for the early phases (where we do not take into account 
radiative feedback). 

• During the embedded phase mass segregation can de- 
velop which persists after gas has been removed and the clus- 
ter has become less dense. This suggests that the observed 
mass segregation of open clusters such as the Pleiades is a 
remnant rather than the result of later secular evolution. 

• Mass segregation followed by partial dissolution of the 
cluster tends to lead to an excess of heavy stars in the clus- 
ter. This is due to the fact that low velocity orbits are pref- 
erentially retained and the energy equipartition that forms 
the basis of mass segregation means that heavier stars pref- 
erentially populate the low velocity orbits. 

• The dependence of the observed distribution of the ra- 
tio of the age of an association and its crossing time (11) is 
consistent with the "resetting of the 11 clock" we find in our 
models. The short crossing times of young dense clusters in- 
crease due to stirring effect of mass loss. Unbound clusters 
are seen to have low fl values. These two effects together 
mean that clusters with both 11 < 1 and 11 > 1 are ex- 
pected to be present as long as the gas dispersal and mass 
loss of young stars happens (before 30 Myr). This is consis- 
tent with what is observed (although more precise modeling 
should determine whether the observed distributions can be 
matched). 

The fact that we find mass segregation and the prefer- 
ential retention of high mass stars has consequences for the 
interpretation of cluster observations. In fact, observed clus- 
ters like the Pleiades show mass segregation. This mass seg- 
regation is usually attributed to the secular evolution in the 
pure N-body regime after the short lived embedded phase. 
For the Pleiades primordi al mass segregation is necessary 
IConverse fc Stahlej 1200^ ). Our models are able to repro- 
duce the expected amount of initial mass segregation for 
the Pleiades cluster (although our models are not specifi- 
cally build to reproduce the Pleaides). The fact that high 
mass stars are better able to withstand the disruptive effect 
of the loss of the natal cloud may lead to an excess number 
of unbound low mass stars in the vicinity of young clusters, 
which may be testable by expanding the search for cluster 
members to include kinematically unbound stars of the right 
age. We have demonstrated that it is important to consider 
both the embedded phase evolution and the gas clearing 
phase of cluster evolution together. Future studies will aim 
to include radiative feedback effects and should also take 
into account more realistic initial conditions for the cluster 
and gas distribution. 
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